library(quantreg)
library(affy)
 data_rtpcr<-read.table("taqman_residual_graph_input.txt",sep="\t",head=T)
data1<-read.table("frma_normalized_taqman1.txt",sep="\t",head=T)
data2<-read.table("frma_normalized_taqman2.txt",sep="\t",head=T)
data3<-read.table("frma_normalized_taqman3.txt",sep="\t",head=T)
data4<-read.table("frma_normalized_taqman4.txt",sep="\t",head=T)
data5<-read.table("frma_normalized_taqman5.txt",sep="\t",head=T)
data6<-read.table("frma_normalized_taqman6.txt",sep="\t",head=T)
data<-read.table("rma_normalizaed_taqman.txt",sep="\t",head=T)



mat_A<-as.matrix(cbind(data1[,1:5],data2[,1:5],data3[,1:5],data4[,1:5],data5[,1:5],data6[,1:5]))
mat_B<-as.matrix(cbind(data1[,6:10],data2[,6:10],data3[,6:10],data4[,6:10],data5[,6:10],data6[,6:10]))
mat_C<-as.matrix(cbind(data1[,11:15],data2[,11:15],data3[,11:15],data4[,11:15],data5[,11:15],data6[,11:15]))
mat_D<-as.matrix(cbind(data1[,16:20],data2[,16:20],data3[,16:20],data4[,16:20],data5[,16:20],data6[,16:20]))

mat_A1<-as.matrix(cbind(data[,1:5],data[,21:25],data[,41:45],data[,61:65],data[,81:85],data[,101:105]))
mat_B1<-as.matrix(cbind(data[,6:10],data[,26:30],data[,46:50],data[,66:70],data[,86:90],data[,106:110]))
mat_C1<-as.matrix(cbind(data[,11:15],data[,31:35],data[,51:55],data[,71:75],data[,91:95],data[,111:115]))
mat_D1<-as.matrix(cbind(data[,16:20],data[,36:40],data[,56:60],data[,76:80],data[,96:100],data[,116:120]))

robustmean <- function(x) { tukey.biweight(x, c = 5, epsilon = 1e-04)}

A<-apply(mat_A,1,robustmean)
B<-apply(mat_B,1,robustmean)
C<-apply(mat_C,1,robustmean)
D<-apply(mat_D,1,robustmean)

A1<-apply(mat_A1,1,robustmean)
B1<-apply(mat_B1,1,robustmean)
C1<-apply(mat_C1,1,robustmean)
D1<-apply(mat_D1,1,robustmean)

model_A<-rq(A~data_rtpcr$A)
pred_A<-predict(model_A,as.data.frame(data_rtpcr$A))
#res_A<-A-pred_A
dval<-mat_A[,1]-pred_A
for(i in 2:ncol(mat_A)){
dval<- cbind(dval,mat_A[,i]-pred_A)
}
res_A<-apply(dval,1,robustmean)



model_B<-rq(B~data_rtpcr$B)
pred_B<-predict(model_B,as.data.frame(data_rtpcr$B))
#res_B<- B-pred_B
dval<-mat_B[,1]-pred_B
for(i in 2:ncol(mat_B)){
dval<- cbind(dval,mat_B[,i]-pred_B)
}
res_B<-apply(dval,1,robustmean)


model_C<-rq(C~data_rtpcr$C)
pred_C<-predict(model_C,as.data.frame(data_rtpcr$C))
#res_C<- C-pred_C
dval<-mat_C[,1]-pred_C
for(i in 2:ncol(mat_C)){
dval<- cbind(dval,mat_C[,i]-pred_C)
}
res_C<-apply(dval,1,robustmean)




model_D<-rq(D~data_rtpcr$D)
pred_D<-predict(model_D,as.data.frame(data_rtpcr$D))
#res_D<- D-pred_D
dval<-mat_D[,1]-pred_D
for(i in 2:ncol(mat_D)){
dval<- cbind(dval,mat_D[,i]-pred_D)
}
res_D<-apply(dval,1,robustmean)


model_A1<-rq(A1~data_rtpcr$A)
pred_A1<-predict(model_A1,as.data.frame(data_rtpcr$A))
#res_A1<-A1-pred_A1
dval<-mat_A1[,1]-pred_A1
for(i in 2:ncol(mat_A1)){
dval<- cbind(dval,mat_A1[,i]-pred_A1)
}
res_A1<-apply(dval,1,robustmean)


model_B1<-rq(B1~data_rtpcr$B)
pred_B1<-predict(model_B1,as.data.frame(data_rtpcr$B))
#res_B1<- B1-pred_B1
dval<-mat_B1[,1]-pred_B1
for(i in 2:ncol(mat_B1)){
dval<- cbind(dval,mat_B1[,i]-pred_B1)
}
res_B1<-apply(dval,1,robustmean)


model_C1<-rq(C1~data_rtpcr$C)
pred_C1<-predict(model_C1,as.data.frame(data_rtpcr$C))
#res_C1<- C1-pred_C1
dval<-mat_C1[,1]-pred_C1
for(i in 2:ncol(mat_C1)){
dval<- cbind(dval,mat_C1[,i]-pred_C1)
}
res_C1<-apply(dval,1,robustmean)


model_D1<-rq(D1~data_rtpcr$D)
pred_D1<-predict(model_D1,as.data.frame(data_rtpcr$D))
#res_D1<- D1-pred_D1
dval<-mat_D1[,1]-pred_D1
for(i in 2:ncol(mat_D1)){
dval<- cbind(dval,mat_D1[,i]-pred_D1)
}
res_D1<-apply(dval,1,robustmean)


library("geneplotter")
require("RColorBrewer")

 pdf("scatterplot_frma_vs_rma_sample_A.pdf")
plot(log(data_rtpcr$A),res_A,  pch=".",xlab="RTPCR", ylab="RESIDUALS",col="red")
lines(lowess(log(data_rtpcr$A),res_A),col="red")
points(log(data_rtpcr$A),res_A1,  pch='.',col="green")
lines(lowess(log(data_rtpcr$A),res_A1), col="green")
title("FRMA vs RMA Scatter Plot:Sample A, 30 replicates")
dev.off()

pdf("scatterplot_frma_vs_rma_sample_B.pdf")
plot(log(data_rtpcr$B),res_B,  pch=".",xlab="RTPCR", ylab="RESIDUALS",col="red")
lines(lowess(log(data_rtpcr$B),res_B),col="red")
points(log(data_rtpcr$B),res_B1,  pch='.',col="green")
lines(lowess(log(data_rtpcr$B),res_B1), col="green")
title("FRMA vs RMA Scatter Plot:Sample B, 30 replicates")
dev.off()

pdf("scatterplot_frma_vs_rma_sample_C.pdf")
plot(log(data_rtpcr$C),res_C,  pch=".",xlab="RTPCR", ylab="RESIDUALS",col="red")
lines(lowess(log(data_rtpcr$C),res_C),col="red")
points(log(data_rtpcr$C),res_C1,  pch='.',col="green")
lines(lowess(log(data_rtpcr$C),res_C1), col="green")
title("FRMA vs RMA Scatter Plot:Sample C, 30 replicates")
dev.off()

pdf("scatterplot_frma_vs_rma_sample_D.pdf")
plot(log(data_rtpcr$D),res_D,  pch=".",xlab="RTPCR", ylab="RESIDUALS",col="red")
lines(lowess(log(data_rtpcr$D),res_D),col="red")
points(log(data_rtpcr$D),res_D1,  pch='.',col="green")
lines(lowess(log(data_rtpcr$D),res_D1), col="green")
title("FRMA vs RMA Scatter Plot:Sample D, 30 replicates")
dev.off()

pdf("Histogram_D.pdf")
h1 <- hist(res_D,breaks=seq(-35,10,0.1))
h2 <- hist(res_D1,breaks=seq(-35,10,0.1))
hh <- cbind(h1$density,h2$density)
colnames(hh) <- c("fRMA","RMA")
barplot2(t(hh),beside=TRUE,col=c("red","blue"),legend=c("RMA","fRMA"),main="RMA vs fRMA residuals",ylab="Counts",xlab="residual Value",xlim=c(800,1200),plot.grid=TRUE)
title("FRMA vs RMA Histogram Plot:Sample D, 30 replicates")
dev.off()

pdf("Histogram_A.pdf")
h1 <- hist(res_A,breaks=seq(-20,10,0.1))
h2 <- hist(res_A1,breaks=seq(-20,10,0.1))
hh <- cbind(h1$density,h2$density)
colnames(hh) <- c("fRMA","RMA")
barplot2(t(hh),beside=TRUE,col=c("red","blue"),legend=c("RMA","fRMA"),main="RMA vs fRMA residuals",ylab="Counts",xlab="residual Value",xlim=c(200,800),plot.grid=TRUE)
title("FRMA vs RMA Histogram Plot:Sample D, 30 replicates")
dev.off()

pdf("Histogram_B.pdf")
h1 <- hist(res_B,breaks=seq(-40,10,0.1))
h2 <- hist(res_B1,breaks=seq(-40,10,0.1))
hh <- cbind(h1$density,h2$density)
colnames(hh) <- c("fRMA","RMA")
barplot2(t(hh),beside=TRUE,col=c("red","blue"),legend=c("RMA","fRMA"),main="RMA vs fRMA residuals",ylab="Counts",xlab="residual Value",xlim=c(950,1350),plot.grid=TRUE)
title("FRMA vs RMA Histogram Plot:Sample D, 30 replicates")
dev.off()


pdf("Histogram_C.pdf")
h1 <- hist(res_C,breaks=seq(-20,10,0.1))
h2 <- hist(res_C1,breaks=seq(-20,10,0.1))
hh <- cbind(h1$density,h2$density)
colnames(hh) <- c("fRMA","RMA")
barplot2(t(hh),beside=TRUE,col=c("red","blue"),legend=c("RMA","fRMA"),main="RMA vs fRMA residuals",ylab="Counts",xlab="residual Value",xlim=c(400,800),plot.grid=TRUE)
title("FRMA vs RMA Histogram Plot:Sample D, 30 replicates")
dev.off()








